{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Expand With Interpolators <a href=\"https://mybinder.org/v2/gh/InsightSoftwareConsortium/SimpleITK-Notebooks/master?filepath=Python%2F20_Expand_With_Interpolators.ipynb\"><img style=\"float: right;\" src=\"https://mybinder.org/badge_logo.svg\"></a>\n",
    "\n",
    "This notebook demonstrates the different interpolators available in SimpleITK available for image resampling.  Their effect is demonstrated on the <a href=\"http://www.cs.cornell.edu/~srm/publications/Vis94-filters-abstract.html\">Marschner-Lobb</a> image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import SimpleITK as sitk\n",
    "import numpy as np\n",
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def myshow(img, title=None, margin=0.05):\n",
    "    if img.GetDimension() == 3:\n",
    "        img = sitk.Tile(\n",
    "            (\n",
    "                img[img.GetSize()[0] // 2, :, :],\n",
    "                img[:, img.GetSize()[1] // 2, :],\n",
    "                img[:, :, img.GetSize()[2] // 2],\n",
    "            ),\n",
    "            [2, 2],\n",
    "        )\n",
    "\n",
    "    aimg = sitk.GetArrayViewFromImage(img)\n",
    "\n",
    "    xsize, ysize = aimg.shape\n",
    "\n",
    "    dpi = 80\n",
    "\n",
    "    # Make a figure big enough to accommodate an axis of xpixels by ypixels\n",
    "    # as well as the ticklabels, etc...\n",
    "    figsize = (1 + margin) * ysize / dpi, (1 + margin) * xsize / dpi\n",
    "\n",
    "    fig = plt.figure(figsize=figsize, dpi=dpi)\n",
    "    # Make the axis the right size...\n",
    "    ax = fig.add_axes([margin, margin, 1 - 2 * margin, 1 - 2 * margin])\n",
    "\n",
    "    t = ax.imshow(aimg)\n",
    "    if len(aimg.shape) == 2:\n",
    "        t.set_cmap(\"gray\")\n",
    "    if title:\n",
    "        plt.title(title)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def marschner_lobb(size=40, alpha=0.25, f_M=6.0):\n",
    "    img = sitk.PhysicalPointSource(\n",
    "        sitk.sitkVectorFloat32, [size] * 3, [-1] * 3, [2.0 / size] * 3\n",
    "    )\n",
    "    imgx = sitk.VectorIndexSelectionCast(img, 0)\n",
    "    imgy = sitk.VectorIndexSelectionCast(img, 1)\n",
    "    imgz = sitk.VectorIndexSelectionCast(img, 2)\n",
    "    del img\n",
    "    r = sitk.Sqrt(imgx**2 + imgy**2)\n",
    "    del imgx, imgy\n",
    "    pr = sitk.Cos((2.0 * math.pi * f_M) * sitk.Cos((math.pi / 2.0) * r))\n",
    "    return (1.0 - sitk.Sin((math.pi / 2.0) * imgz) + alpha * (1.0 + pr)) / (\n",
    "        2.0 * (1.0 + alpha)\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "myshow(marschner_lobb())\n",
    "myshow(marschner_lobb(100))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = marschner_lobb()\n",
    "ml = ml[:, :, ml.GetSize()[-1] // 2]\n",
    "myshow(sitk.Expand(ml, [5] * 3, sitk.sitkNearestNeighbor), title=\"nearest neighbor\")\n",
    "myshow(sitk.Expand(ml, [5] * 3, sitk.sitkLinear), title=\"linear\")\n",
    "myshow(sitk.Expand(ml, [5] * 3, sitk.sitkBSpline3), title=\"b-spline order 3\")\n",
    "myshow(sitk.Expand(ml, [5] * 3, sitk.sitkGaussian), title=\"Gaussian\")\n",
    "myshow(\n",
    "    sitk.Expand(ml, [5] * 3, sitk.sitkHammingWindowedSinc),\n",
    "    title=\"Hamming windowed sinc\",\n",
    ")\n",
    "myshow(\n",
    "    sitk.Expand(ml, [5] * 3, sitk.sitkBlackmanWindowedSinc),\n",
    "    title=\"Blackman windowed sinc\",\n",
    ")\n",
    "myshow(\n",
    "    sitk.Expand(ml, [5] * 3, sitk.sitkCosineWindowedSinc), title=\"Cosine windowed sinc\"\n",
    ")\n",
    "myshow(\n",
    "    sitk.Expand(ml, [5] * 3, sitk.sitkWelchWindowedSinc), title=\"Welch windowed sinc\"\n",
    ")\n",
    "myshow(\n",
    "    sitk.Expand(ml, [5] * 3, sitk.sitkLanczosWindowedSinc),\n",
    "    title=\"Lanczos windowed sinc\",\n",
    ")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
